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I | Abstract. We reconstruct the polarization sector of a bright polarization squeezed beam 

("^ starting from a complete set of Stokes measurements. Given the symmetry that underlies 

O | the polarization structure of quantum fields, we use the unique SU(2) Wigner distribution 

I to represent states. In the limit of localized and bright states, the Wigner function can be 

_i approximated by an inverse three-dimensional Radon transform. We compare this direct 

Jlj reconstruction with the results of a maximum likelihood estimation, finding an excellent 

_. agreement. 
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1. Introduction 



Polarization of light is a robust characteristic that can be efficiently manipulated using modest 

equipment without introducing more than marginal losses. It is thus not surprising that this 
S-J is often the variable of choice to encode quantum information, as one can convince oneself 

^ by looking at some recent cutting-edge experiments, including quantum key distribution JT], 

y—{ quantum dense coding |2|, quantum teleportation |3|, rotationally invariant states [4|, phase 

super-resolution [5|, and weak measurements (6). 

In the discrete-variable regime of single, or few, photons, one is mostly interested into 

two-mode states, which for all practical purposes can be regarded as a spin system [7 , 8|. 

As a result, the polarization state can be determined from correlation functions of different 

orders |[9]jT9). Given the small dimensionality of the Hilbert space involved, the state 

reconstruction can be readily performed. 

In the continuous-variable case, polarization properties are exploited for an expedient 



generation, manipulation, and measurement of nonclassical light. Polarization squeezing [20- 
23 1, which has been observed in numerous experiments |24-28|, is perhaps the most 
tantalizing illustration. Full Stokes polarimetry |29] is the method employed by the majority 
of the practitioners in this area. 

However, the reconstruction in this limit is a touchy business and requires special care. 
The origin of the problem can be traced back to the fact that the characterization of the 
polarization state by the whole density operator is superfluous, because it contains much 
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more than polarization information. This redundancy can easily be handled for low number 
of photons, but becomes a real hurdle for highly excited states. An adequate solution has 
been proposed recently: it suffices with a subset of the density matrix that has been called the 
"polarization sector" |30|31| or the "polarization density operator" (32} . Its knowledge allows 



for a complete specification of the state on the Poincare sphere (actually on a set of nested 
spheres that can be appropriately called the Poincare space). The technique was devised by 
Karassiov and coworkers 1 33 -35) and implemented experimentally in [36]. 



In this paper we present a comprehensive treatment of polarization tomography. As with 
any reliable quantum tomographical scheme, we need to supply three key ingredients |37|: 
the availability of a tomographically complete measurement, a suitable representation of the 
quantum states, and a robust algorithm for inverting the experimental data. In this respect, we 
use a standard Stokes scheme that implements the first ingredient in a very simple way; for the 



second, we resort to the well-known SU(2) Wigner distribution [38-45], and finally, we prove 
that the inversion of the data in terms of that Wigner function is an inverse three-dimensional 
(3D) Radon transform. 

To support the experimental feasibility of our scheme, we carry out the full tomography 
of a bright polarization-squeezed state generated in a Kerr medium [46]. The reconstruction 
is accomplished in three-different ways: by the direct inversion of the Radon transform, by a 
maximum-likelihood estimation and, finally, by a Gaussian approximation. The final results 
are compared and the sources of uncertainty are analyzed. 

2. Setting the theoretical foundations 

2.1. Polarization structure of quantum fields 

A satisfactory description of the polarization structure of quantum fields and the 
corresponding observables that specify this structure is of paramount importance for our 
purposes. 

We restrict our attention to the case of a monochromatic plane wave (the formalism 
can be easily extended to more involved multimode wavefronts [47 , 48 1), which we assume to 
propagate in the z direction, so its electric field lies in the xy plane. Under these conditions, we 
are dealing with a two-mode field that can be fully characterized by two complex amplitude 
operators. They are denoted by &h and ay, where the subscripts H and V indicate horizontally 
and vertically polarized modes, respectively. The commutation relations of these operators are 

[a k ,4] = S U , k,ee{H,V}. (2.1) 

The description is greatly simplified if we use the Schwinger representation J49][50) 

J\ = \{a H a v + dyd H ) , f 2 = | (feay - cL^ay) , h = \{d* H a H - a v a v ) , (2.2) 

together with the total number 

N = a H an + dydy . (2.3) 

These operators coincide, up to a factor 1/2, with the Stokes operators [51], whose average 
values are precisely the classical Stokes parameters (52J. Using equation ( |2.1| >, one 
immediately notices that J = {J\,J2 1 J-i) satisfy the commutation relations of the su(2) algebra 

[Jk,Je] = iZMmJm , (2.4) 

where ejdm is the Levi-Civita fully antisymmetric tensor. This noncommutability precludes 
the simultaneous exact measurement of the physical quantities they represent. Among 
other consequences, this implies that no field state (apart from the vacuum) can have 
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sharp nonfluctuating values of all the operators J simultaneously. This is expressed by the 
uncertainty relation 

A 2 J = A 2 j\ + A 2 / 2 + A 2 / 3 > (N) /2 , (2.5) 

where the variances are given by A 2 ft = (jf) — (fj) 2 . In other words, the electric vector of a 
monochromatic quantum field never traces out a definite ellipse. 

In classical optics, the total intensity is a well-defined quantity and the Poincare sphere 
appears thus as a smooth surface with radius equal to that intensity. In quantum optics we 
have 

^+-/!+/!=(!)(f + l), (2-6) 

and, as fluctuations in the number of photons are unavoidable (leaving aside photon-number 
states), we are forced to talk of a three-dimensional Poincare space (with axis J\, Ji and Jj,) 
that can be envisioned as foliated in a set of nested spheres with radii proportional to the 
different photon numbers that contribute significantly to the state. 

The Hilbert space Jf of these two-mode fields has a convenient orthonormal basis in the 
form of Fock states for both polarization modes, namely |n#,ny), However, since 

[N,J]=0, (2.7) 

each subspace with a fixed number of photons N must be handled separately. In other words, 
in the previous onion-like picture of the Poincare space, each shell has to be addressed 
independently. This can be emphasized if instead of the Fock basis, we employ the relabeling 

\J,m) = \nj] =J + m,riy = J — m) , (2.8) 



According to (2.6 1, we have that J = N/2 and this basis can be also seen as the common 
eigenstates of {J Z ,J^}. In this way, for each fixed J (i.e., fixed number of photons AO, m runs 
from —J to / and these states span a (27+ 1) -dimensional subspace wherein J acts in the 
usual way (in units h = 1) 

f±\J,m) = y/j(J+l)-m(m±l)\J,m±l), 

(2.9) 

jT,\J,m) = m\J,m) , 

with./± =Ji±J2- 

It is clear from all this previous discussion that the moments of any energy-preserving 
observable (such as J) do not depend on the coherences between different subspaces. The only 
accessible information from any state described by the density matrix p is thus its polarization 
sector, which is specified by the block-diagonal form 

P P oi = 0P (7) (2-10) 

j 

where p^ J > is the reduced density matrix in the J subspace. Any p and its associated block- 
diagonal form p po i cannot be distinguished in polarization measurements (and, accordingly, 
we drop henceforth the subscript pol). This is consistent with the fact that polarization and 
intensity are, in principle, independent concepts: in classical optics the form of the ellipse 
described by the electric field (polarization) does not depend on its size (intensity). 
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2.2. Polarization squeezing and the dark plane 



The variances of the angular-momentum operators (2.2 1 are not independent, for they are 
constrained by 

A 2 J k A 2 f e >e khn \(J m )\ 2 . (2.11) 

It is always possible to find pairs of maximally conjugate operators for this uncertainty 
relation. This is equivalent to establishing a basis in which only one of the operators ( |2.2| i has 
a nonzero expectation value, say (J k ) = (J() — and (J m ) ^ 0. The only nontrivial Heisenberg 
inequality reads thus 

A 2 / A .A 2 ^>|(4)| 2 . (2.12) 



Polarization squeezing can then be sensibly defined as 1 20 - 23 1: 

A 2 A.<|(7;,)|<A 2 ^. (2.13) 

The choice of the conjugate operators {J k ,J(} is by not means unique: there exists an infinite 
set {J^(9),J^(9 + 7C/2)} that are perpendicular to the state classical excitation 7 m , for which 
(J ±(9)) = for all 9. All these pairs exist in the J k -Je plane, which is called the "dark 
plane" because it is the plane of zero mean intensity. We can express a gen eric J± (9) as 



J ±(9) — J k cos 9 +Ji sin0, 9 being an angle defined relative to J k . Condition (2.13 1 is then 
equivalent to 

A 2 / ± (0 sq ) < i|(JV)| < A 2 / ± (e sq + 7T/2), (2.14) 

where ,/x(0sq) is the squeezed parameter and J±(9 S q + 7t/2) the antisqueezed parameter. 

In the experiments presented in this paper, a focal role will be played by the example 
in which the horizontal and vertical modes have the same amplitude but are phase shifted by 
k/2: (an) = i{&y) = ia/y/2, a being a real number. This light is circularly polarized and 
fulfills (J i ) = (fj) = 0, (J2) —a 2 . It is advantageous to work in the circular polarization basis, 
whose right (+) and left (— ) amplitudes are given in terms of the linear ones by 

a± — —=(a H ±ia v ) . (2.15) 

v2 

In this manner, (d + ) = a and (a_) = 0. The operators in the J1-J1 plane correspond to the 
quadrature operators of the dark left-polarized mode. In fact, expressing the fluctuations of J 
in terms of the noise of the circularly polarized modes 8d± and assuming |(5a±)| CHwe 
find (53) 

8J\(9) = a8X-(9) = a[8X H {9) + 8X V {9 +n/2)\ , (2.16) 

where X[ = (a,e~' e +d j e' e )/\/2 are the rotated quadratures for the /th amplitude. On the 
other hand, we have that 

8N = a{8a+ + 8dl) = a8X+, (2.17) 

and the intensity exhibits no dependence on the dark mode. In consequence, the condition 
( |2.14[ ) can be recast for this example as 

A 2 7 L (0)<|(a)| 2 & A 2 X_{9)<1, (2.18) 

that is, polarization squeezing is equivalent to vacuum squeezing in the orthogonal 
polarization mode. 

In the dark-plane measurements, the beam is divided equally between two 
photodetectors. Such measurements are then identical to balanced homodyne detection: 
the classical excitation is a local oscillator for the orthogonally polarized dark mode. The 
phase between these modes is varied by rotating the measurement through the dark plane, 
allowing a full characterization of the noise properties. This is a unique feature of polarization 
measurements and has been used in many experiments [54-58 1. 
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2.3. Polarization Wigner function 

The structure discussed so far highlights that SU(2) is the symmetry group for polarization. To 
provide an appropriate phase-space description of states, we take advantage of the pioneering 
work of Stratonovich [38 1 and Berezin {39) , who introduced quasi-probability distribution 
functions on the sphere satisfying all the proper requirements. This construction was later 
generalized by others [40-45] and has proved to be very useful in visualizing properties of 
spinlike systems | |59] - |63| 

To gain physical insights into this approach, let us start by representing the density matrix 
with respect to the polarization basis. Instead of using directly the states {\J,m)}, it is more 
convenient to write such an expansion as 

P {J) =t ip$t$, (2-19) 

K=0q=-K 



where the irreducible tensor operators T% are [64 



1k " "V 27- 



j £ CJi Kq \J,m')(J,m\, (2-20) 

m,m'=—J 



with CJ™ K being the Clebsch-Gordan coefficients that couple a spin 7 and a spin K 
(0 < K < 27) to a total spin 7. These tensor operators have the right transformation properties 
under rotations and they indeed constitute the most suitable orthonormal basis 

nT$f^} = 8 J j,8 KK ,8 qq ,. (2.21) 

Although at first sight they might look a little bit intimidating, they are nothing but the 
multipoles used in atomic physics |65|: one can check that 

■>(■*) _ 1 jj f(J) = I 3 

V27TT •* V(2J+l)(J + l)J' 



f oc \' = ^r^ 1 f u=\ h^ n/„ u r J i « = ± > z - ( 2 - 22 ) 



and similarly the 7^ can be related to the Kth power of the generators (|2.2|i. Accordingly, 



the expansion coefficients 

pg=Tr[pWr]f] (2.23) 

are known as state multipoles. 

The Wigner function associated with the state ( |2.19 i is 

W {J) (9,<i>)=Tr[p {J) A {J) (9,<l))], (2.24) 

where A^> ( 9 , ) is the Wigner kernel 



TJ K 



A(y) ( ^) = \/y7ZTl I Y^(e,t)t$, (2.25) 

' /J + 1 K=0q=-K 

and Ykc/(9,(J)) are the spherical harmonics. This kernel is unitary and satisfies the 
normalization conditions 

Tr[AW(0,A)] = l, 2 £±1 f dO.W\6,6) = l. (2.26) 

An jyi 

The integral extends over the unit sphere ^ 2 and dQ. is the invariant measure therein, namely, 
dQ. = sin9d9d(j). 
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From ( |2.24| > and the properties of the irreducible tensors, one can immediately express 
the Wigner function in the very suggestive form 

w w (M)=£ t pSW*)' (2 - 27) 

K=Qcj=-K 

which clearly shows that determining this Wigner function is tantamount to the knowledge of 
all the state multipoles. (i.e., all the moments of the Stokes parameters). 

One can obtain the marginal of W^ 0,0) once summed over all the values of J 

m0)=I^pW«(#,<p), (2.28) 

where the factor has been introduced to ensure the proper normalization. 

In the above-mentioned example of a strong circularly polarized state, we can consider 
that the sphere can locally be replaced by its tangent plane since J ~ a 2 . Using simple 
geometrical relations between the coordinates (0,0) and the Cartesian coordinates (q,p) in 
that tangent plane, we get 

W(d,<l>)~aW(q,p), (2.29) 

which confirms that the dark plane is equivalent to the standard phase space for continuous 
variables. 

In the limit of large photon numbers the representation ( |2.24[ ) is not very useful. In such 
a case, a remarkably effective approximation is given by (66) 

AW(0,0)~(-l)' / exp(-^n-J), (2.30) 

where n = (cos sin , sin sin , cos ) is the unitary vector in the direction (0,0). 

2.4. Tomograms and tomographic inversion 

A general polarimetric apparatus consists of a half-wave plate, with axis at angle a, followed 
by a quarter- wave plate at angle p\ For fixed values of the angles (<X,j3) of the wave plates, 
the selected direction in the Stokes space is 

= ?r/2-2j3, 0=2j3-4a. (2.31) 

The polarization transformations performed by the wave plates can be represented by Ji, 
which generates rotations about the direction of propagation, and J$, which generates phase 
shifts between the modes. Their joint action is given by the operator 

D(0^)=e iej ' 2 e^ J '\ (2.32) 

which describes displacements over the sphere. After that, a polarizing beam splitter projects 
onto the basis \J,m). 

In physical terms, the wave plates transform the input polarization allowing the 
measurement of different Stokes parameters by the projection onto the basis \J,m). This 
can be modeled by 

t$ = \J,m)(J,m\, (2-33) 

so that w m = Tr[pn,„ ] is the probability of detecting n# = J + m photons in the horizontal 
mode and simultaneously ny = J — m photons in the vertical one. Of course, when the total 
number of photons is not measured and only the difference m is observed, it reduces to 

n,„= £ \J,m){J,m\. (2.34) 

J=\m\ 
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The experimental histograms recorded for each direction (6,<p) correspond to the 
tomographic probabilities 

(y) (0,0)=Tr[pn£ /) (0^)]=Tr[pD(0,0)nL y) D t (0,0)]. (2.35) 



ii', 



JT{x) = ^ r j d©sm 2 ( k yJ«- M ». (2.39) 



The reconstruction in each (27+ 1) -dimensional invariant subspace can now be carried out 
exactly since it is essentially equivalent to a spin J |67] - [70) . One can proceed in a variety of 
ways, but perhaps the simplest one is to look for an integral representation of the tomograms 
( |2.35[ ); as soon as one realizes that 

n;,, J (0,0) = — / rffflexp totf-n-ro) , (2.36) 

2k Jo 

the tomograms read as 

w£P(e,0) = -5- / K dcoTv[p iJ] exp(/0)j-n)] e - ! ' fflm , (2.37) 

27C Jo 

that is, they appear as the Fourier transform of the characteristic function for the observable 

J • n. After some direct manipulations, we find that 

P {J) = J- L / ^n'wL 7) (n')jr(J-n'~m), (2.38) 

where dn! indicates integration over the unit sphere and the kernel ,$f(x) is 

2/+1 [ ln , . 2 ( co- 
da} sir ' 
/o 

Although ( 2.38) is a formal solution, it is handier to map this density matrix onto the 

corresponding Wigner function, for which we need to compute Tr[J^(J • n' — m) A ^ (0,0)]. 

When J is large enough, we can replace the Wigner kernel by its approximate expression 

( |230| >, getting 

Tr[jr(J-n'-m)A( J )(0,</.)] = {-\) J ^^ J da} sin 2 (|) ^'""^(a)') , (2.40) 

where ^y(fi)') is the character of the (27 + 1) -dimensional representation of the SU(2) 
group [64] and (O 1 is given by 

cos (y) =n ' n ' sin (f) • (2 - 41) 

For J ^> 1, m can be taken as a continuous variable. Replacing the sum by an integral, 
integrating by parts and taking into account that for localized states n-n' ~ 1, the Wigner 
function simplifies to 

W(J,6,<1>) = 2 ^ f dmf dn' d -^ n ^8(m~Jn-n'), (2.42) 

4k z J-oo Jy 2 dm z 

where we have included J as an argument to stress that it must be treated as continuous. The 
reconstruction in this limit turns out to be equivalent to an inverse Radon transform of the 
measured tomograms. 
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3. Experiment 

3.1. The setup 

To validate our approach, we perform the tomography of a polarization squeezed state, 
generated in a polarization-maintaining optical fibre through the nonlinear Kerr effect |46|. 
The setup is shown in figure 1 . The light source is a shot-noise limited ORIGAMI laser from 
Onefive GmbH emitting 220 fs pulses at a repetition rate of 80 MHz and centered at 1560 nm. 
The light is fed into a 13 m-long polarization-maintaining birefringent fibre (3M FS-PM- 
7811, 5.6 jiim mode-field diameter) so that quadrature squeezed states are simultaneously and 
independently generated in both polarization modes. The strong birefringence of the fibre 
(beat length 1.67 mm) causes a delay between the emerging quadrature-squeezed pulses. We 
precompensate for this delay in an unbalanced Michelson-like interferometer placed before 
the fibre. A small part (0.1 %) of the fibre output serves as the input to a control loop to 
maintain the relative phase between the exiting pulses locked to Tt/2, so the light is circularly 
polarized. 

The quantum state is detected with a Stokes measurement, as sketched in the previous 
section. The two detectors (InGaAS PIN photodiodes, custom-made by Laser Components 
GmbH with 98 % quantum efficiency at DC) are balanced and have a sub-shot noise resolution 
at a frequency range between 5 MHz and 30 MHz. Each detector has two separate outputs: 
DC, providing the average values of the photocurrents, and AC, providing the photocurrents 
amplified in radio-frequency (RF) spectral range. The RF currents of the photodetectors are 
mixed with an electronic local oscillator at 12 MHz, amplified (FEMTO DHPVA-100), and 
digitized by an analog/digital exit converter (Gage CompuScope 1610) at 10 Megasamples 
per second with a 16-bit resolution and 10 times oversampling. 

The measurements are performed at a pulse energy of 93 pJ. In the dark plane a total 
squeezing of about 3.8 dB is observed. In the orthogonal quadrature, the noise was enhanced 
by several tens of dB due to Guided Acoustic Wave Brillouin Scattering (GAWBS) |7T[{73") . 
In the direction of the classical excitation, the state is expected to be shot-noise limited, since 
the Kerr effect only influences the phase and does not contribute to the photon number. 

To perform the reconstruction, histograms of the Stokes variables are recorded for 
different angles (0,0). This is done by rotating the wave plates with motorized stages 
(OWIS DMT 40-D20-HSM) and scanning one eighth of the Poincare sphere in 8100 steps, 
a measurement that took over eight hours. The unmeasured octants were deduced from 
symmetry. For each setting of the wave plates, the photocurrent noise of both detectors was 



Birefringence compensator 



Stokes measurement 




Figure 1. Setup for efficient polarization squeezing generation and the conesponding Stokes 
measurement apparatus. 
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Figure 2. Measured histograms of the difference current of the two detectors for various 
measurement directions on the Poincare sphere.. Note the different scales on both plots. 
Histograms 1, 2, and 3 are in the dark plane, while histogram 4 is at the classical mean value. 
The histograms corresponding to electronic and shot noise are also shown. 



simultaneously sampled 0.5 x 10 6 times. Noise statistics of the detectors difference current 
were acquired in histograms with 751 bins. Additionally, the optical intensity was recorded. 

In figure 2 we show typical histograms at different angles on the Poincare sphere. As 
the widths largely vary from squeezing to antisqueezing ranges, there are two plots in which 
the amplitude scale differs in more than one order of magnitude. The histograms labeled 1, 2 
and 3 are measured in the dark plane. Tomogram 1 denotes the angle of maximum squeezing, 
while 3 corresponds to the antisqueezing. Tomogram 4 is at the classical mean value, where 
the measured noise is almost shot-noise limited. Due to the high number of samples, the 
measured histograms are smooth and, at the same time, the number of bins makes it possible 
to resolve the large dynamical range of amplitudes, so no data interpolation was needed. We 
also plot histograms showing the electronic noise and the shot noise. 

For all these histograms we have performed a Gaussianity check, using the Kolmogorov- 
Smirnov and the % 2 tests, as well as the Kullback-Leibler divergence [74]. We can conclude 



that all the histogramas are Gaussian within confidence levels ranging from 95 % to 98 %. 



3.2. Experimental reconstruction 



As clearly expressed in (2.42 1, for high photon numbers the tomography turns out to be 



equivalent to an inverse Radon transform of the measured histograms. In practice, this one- 
step 3D reconstruction is very demanding in computational resources. Therefore, we divide 
the process into two steps: first, a set of 2D projections is reconstructed from the recorded 
histograms; subsequently, the Wigner function is slice-wise generated from those projections 
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Figure 3. (Right) Isocontour surface of the level 1/e (from the maximum) of the Wigner 
function W (J, 9 , ) for the polarization squeezed state generated in the setup of figure 1 . (Left) 
Sections of the Wigner function through the three coordinate planes. In blue we show the 
isotropic section for a coherent state, which we use as unit for all the plots. 



(to which we apply a Hamming filter to smooth the noise). The symmetry of the state is used 
as a prior information to reduce the range of measured angles to an octant. This minimizes 
the systematic errors stemming from imperfections of the polarization optics. 

In figure 3 (right panel) we show the final result of the inverse Radon transform for our 
polarization squeezed state. More concretely, we plot an isocontour surface of W(J, 9,<j>) = 
constant (with the constant being 1/e from the maximum) in the Poincare space having J\, 
J-i, and Jj as orthogonal axes. As coordinate units we use the shot noise set by a coherent 
state. The ellipsoidal shape of the state is clearly visible. The center of the ellipsoid is far 
away from the origin, since we have 10 x 10 11 photons per measurement time (using 1.9 MHz 
resolution bandwidth). The antisqueezed direction of the ellipsoid is dominated by excess 
noise stemming largely from GAWBS, as we have already mentioned. 

In the left panel of figure 3 we sketch density plots of the projections on the coordinate 
planes of the previous Wigner function (including the particular case of a coherent state). The 
contours agree with the 3.8 ± 0.3 dB squeezing that was directly measured from the variances. 
The projections on the planes J1-J2 and 72-^3 show an additional spreading of the state in the 
J2 direction caused by the imperfect polarization contrast in the measurement setup that mixes 
some of the antisqueezing on the J2 direction. 

This Radon reconstruction requires a large set of measured data to get a reasonably 
accurate representation of the state. There are two main reasons for this: integrals are 



approximated by finite sums (in our case, we used 751 bins in 91 steps) and the kernel (2.39 1 
is singular, so some ad hoc filtering of the raw data is needed. Acquiring such large data sets 
may be unwise, for it demands long measurement times. Ensuring the proper stability of the 
setup is thus essential and might be difficult depending on the quantum state measured. 

This limitation may be circumvented by adopting a statistically-motivated method, such 
as the maximum likelihood (ML) [75]. In our case, the relation between the Wigner function 
W and the tomograms w can be written as a system of linear equations 



W j 



I> ; ^, 



(3.1) 



where the subscripts in Wk and w i is a shorthand notation for the respectives coordinates. 
The coefficients Cjk can be interpreted as the overlap of the jth projector with the kth 



volume element of the Wigner function and can be readily determined from equations (2.27 1 



and (2.35 1. The most likely Wigner function is then found by minimizing the Kullback- 



Leibler divergence between the normalized vectors of the computed tomograms wj and 
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Figure 4. Sections of the isocontour of the Wigner function in figure 3 through the coordinate 
axes, for the three different reconstruction techniques used. In grey, direct Radon transform, 
in green ML method with nine settings for the angles of the wave plates and in dotted lines the 
results of a Gaussian ML approximation. 



recorded ones wj. Technically, this can be achieved by the iterative expectation-maximization 



algorithm [76-78 1 



(LjWjjiLjCjkjywj 

which converges monotonously to the ML estimate from any strictly positive initial vector 

The significantly greater stability of the statistical inversion allows us to get 
reconstructions of the same quality but from far smaller data sets. This is illustrated in 
figure 4, where we draw a comparison between Radon and ML methods, although in the latter 
case using only nine different settings of the angles (0,<j>), which amounts to reducing the 
measurements by two orders of magnitude. In other words, the measuring time is shortened 
from eight hours to less than five minutes! This result indicates that the experimental 
characterization of considerably more complicated quantum states with less symmetries 
should still be within the reach of the present measurement setup. 



As we have discussed in section 2.2 the dark plane is of special interest. The theory 
shows that the reconstruction therein can be obtained in two different ways: either by 
reconstructing the dark mode directly from the histograms or by calculating projection of 
the 3D Wigner function along the Ji direction. The two results are compared in figure [5] and 
good agreement within the experimental uncertainties is found. Since the Radon transform 
should provide a plausible explanation for all the measured histograms, such a comparison 
may serve as an independent test of the quality of the 3D tomography. 

Finally, the high confidence levels of the Gaussianity tests seems to call for a Gaussian 
ML reconstruction. The Gaussianity is used as a prior information about the signal, which 
helps to reduce drastically the number of free parameters. In this case, the Wigner function is 
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Figure 5. Dark plane reconstructions. Left panel: Reconstruction obtained by integrating the 
Wigner function shown in figure. 3 in the Ji direction. Right panel: ML reconstruction from 
dark-plane histograms. Only nine settings of angles 8 and were used for the ML tomography. 



represented by the 3x3 covariance matrix G: 

W(n)°cexp(--nG- 1 n) , (3.3) 

and the calculated variances Oj = HjGrij are matched (in the ML sense) to the actually 
measured variances [79] . Since G must be positive semidefinite, only six real parameters 
describe the measured system and the problem is highly overdetermined, in consequence, the 
Gaussian state can be obtained from a few histograms. In principle, by comparing Gaussian 
reconstructions based on different subsets of measured data, various imperfections of the 
setup, such as instabilities and biases, can be detected. 
The matrix G turns out to be 

/ 3.0920 xlO 2 -1.1931 -2.0160 \ 

G= -1.1931 4.4485 x 10" 1 -1.2926 x 10~ 2 , (3.4) 

\ -2.0160 -1.2926 x 10~ 2 1.1511 / 

which once diagonalized gives the principal variances 0.43962, 1.13853, and 309.22177 (in 
shot-noise units). This agrees well with the standard and ML reconstructions, as can be also 
appreciated in figure 4. The Gaussian reconstruction was done without assuming a particular 
orientation or symmetry of the state with respect to the Stokes coordinates. The covariance 
matrix suggests that the misalignment of the principal axis is less than 0.5 degrees within the 
measurement errors, in accordance with the definition of angles adopted in the experiment. 

This Gaussian approach allows for a simple estimate of the errors: just take the 
pseudoinversion of the measurement matrix as a linear model and use the standard theory 
of error propagation. The errors to be propagated are actually the errors in the estimated 
variances for each tomogram, which are found from the % 2 distribution. In addition, we can 
assume that the variances of different tomograms are uncorrelated. 

Taking a 97.5 % confidence interval (which corresponds to three standard deviations), 
the principal variances can be written as 

0.440±0.002, 1.139±0.001, 309.2±0.3. (3.5) 

Note that the relative errors in the two larger variances are roughly the same (~ 0. 1 %), while 
for the smallest variance is four times larger. This is a consequence of the experimental setup: 
the smallest variance is directly revealed only in one of the recorded projections used for the 
reconstruction. 
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4. Concluding remarks 

In summary, we have presented a complete programme for the full polarization tomography 
of quantum states. Using the SU(2) Wigner function, we have provided an exact inversion 
formula in terms of the histograms of a standard Stokes measurement and derived a simplified 
version for very localized, high intensity states which turns out to be an inverse Radon 
transform. As a test of the theory, the reconstruction of an intense polarization squeezed 
state has been performed. A ML reconstruction algorithm has also been presented and has 
been compared to the direct method, thereby yielding an excellent agreement. Of course, the 
technique can be readily used for any other polarization state. 
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